options ls=80;
libname dat 'C:\Work Projects\Timeseries_CI\Final texts numbers pgms and outputs';
/*---------------------------;
-----Total Prescription------;
--------------------------- ;*/

data drug1;set dat.red_war2;
     month=_n_;
     intervention=(month>=36);
     monthafintervention=max(0,month-35);
run;
data drug2;set drug1;
     if 39>=month >=36 then delete;
	 if month >= 36 then month=month-4;
	 monthafintervention=max(month-35,0);
run;
proc autoreg data=drug2;
     model war_other5=month intervention monthafintervention/backstep slstay=0.05 nlag=12 method=ml covb;
	 model war_other5=intervention monthafintervention/backstep slstay=0.05 nlag=12 method=ml covb;
run;
*---------------;
*Delta method;
*---------------;
*Absolute and relative change and 95% CI for level change;
%diffrr(drug2,drg_lv,war_other5,intervention monthafintervention,12,1 0 0,1 1 0,0.05);
*Absolute and relative change with 95% CI six months after the alert;
%diffrr(drug2,drg_sixm,war_other5,intervention monthafintervention,12,1 0 0,1 1 6,0.05);
*---------------------;
*Bootstrapping method;
*---------------------;
options nonotes nosyntaxcheck nosymbolgen;
%boots(3263,0,-311.3794,-21.2870,52,
       25993,%str(1),10000,36,12,
       time intervention1 timeafterintervention1);
data meanvar;
     set simulationresults(firstobs=2);
	 lvrel=intervention1/intercept;
	 *slrel=timeafterintervention1/time;
	 sixmabs=intervention1+6*timeafterintervention1;
     sixmrel=(intervention1+6*timeafterintervention1)/(intercept+41*time);
run;
%macro nm(vr);
proc univariate data=meanvar alpha=0.05;
        var &vr;
	    output out=&vr pctlpre=col pctlpts=2.5,50,97.5 mean=me
        ;
run;
proc print;run;
%mend nm;
*Use median is better;
%nm(lvrel);%nm(sixmrel);

/*-----------------------------------;
-----Acetaminophen Prescription------;
------------------------------------;*/
*---------------;
*Delta method;
*---------------;
data drug1;set dat.red_war2;
     month=_n_;
     intervention=(month>=36);
     monthafintervention=max(0,month-35);
run;
data drug2;set drug1;
     if 39>=month >=36 then delete;
	 if month >= 36 then month=month-4;
	 monthafintervention=max(month-35,0);
run;
proc autoreg data=drug2;
     model war_apap=month intervention monthafintervention/backstep slstay=0.05 nlag=12 method=ml covb;
run;
*Absolute and relative change and 95% CI for level change;
%diffrr(drug2,drg_lv,war_apap,month intervention monthafintervention,12,1 0 0 0,1 0 1 0,0.05);
*Absolute and relative change and 95% CI for trend change;
%diffrr(drug2,drg_slope,war_apap,month intervention monthafintervention,12,0 1 0 0,0 1 0 1,0.05);
*Absolute and relative change with 95% CI six months after the alert;
%diffrr(drug2,drg_sixm,war_apap,month intervention monthafintervention,12,1 41 0 0,1 41 1 6,0.05);
*---------------------;
*Bootstrapping method;
*---------------------;
options nonotes nosyntaxcheck nosymbolgen;
%boots(2321,8.2995,-310.9912,-22.1976,52,
       24017,%str(1),10000,36,12,
       time intervention1 timeafterintervention1);
data meanvar;
     set simulationresults(firstobs=2);
	 lvrel=intervention1/intercept;
	 slrel=timeafterintervention1/time;
	 sixmabs=intervention1+6*timeafterintervention1;
     sixmrel=(intervention1+6*timeafterintervention1)/(intercept+41*time);
run;
%macro nm(vr);
proc univariate data=meanvar alpha=0.05;
        var &vr;
	    output out=&vr pctlpre=col pctlpts=2.5,50,97.5 mean=me
        ;
run;
proc print;run;
%mend nm;
*Use median is better;
%nm(lvrel);%nm(slrel);%nm(sixmrel);